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We study structure and fluid-phase behaviour of a binary mixture of hard spheres 
(HSs) and hard spherocylinders (HSCs) in isotropic and nematic states using the 
NPnAT ensemble Monte Carlo (MC) method in which a normal pressure tensor 
component is fixed in a system confined between two hard walls. The method allows 
one to estimate the location of the isotropic-nematic phase transition and to observe 
the asymmetry in the composition between the coexisting phases, with the expected 
increase of the HSC concentration in the nematic phase. This is in stark contrast 
with the previously reported MC simulations where a conventional isotropic NPT 
ensemble was used. We further compare the simulation results with the theoretical 
predictions of two analytic theories that extend the original Parsons-Lee theory using 
the one-fluid and the many-fluid approximation [Malijevsky at al J. Chem. Phys. 
129 , 144504 (2008)]. In the one-fluid version of the theory the properties of the 
mixture are mapped on an effective one-component HS system while in the many- 
fluid theory the components of the mixtures are represented as separate effective HS 
particles. The comparison reveals that both the one- and the many-fluid approaches 
provide a reasonably accurate quantitative description of the mixture including the 
predictions of the isotropic-nematic phase boundary and degree of orientational order 
of the HSC-HS mixtures. 
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I. INTRODUCTION 


Advanced materials formed by the self-assembly of non-spherical building blocks has 
experienced an unprecedented growth due to recent advances in experimental techniques to 
create nano and colloidal particles with almost any imaginable shape.-^- Functional materials 
can be engineered by tailoring the properties of the individual building blocks.-*^ Colloidal 
particles are particularly attractive as building blocks for the design of mesoscale materials, 
which are difficult to fabricate by using chemical synthesis, as their interactions can also 
be modulated by modifying both the surface chemistry of the particles as well and the 
properties of the solvent media.- It is possible to tune the interactions of either steric- 
stabilised or charge-stabilised colloidal particles as nearly as hard body-like by matching 
the index of refraction of both particles and solvent.- Moreover, the self-assembly of these 
systems can also be controlled by the aid of external forces such as magnetic and electric 
fields, gravity^, and even the use of geometrical confinement.-"— We refer to these processes 
in general as directed self-assembly.— 

Anisotropic particles can exhibit many fascinating structures in bulk and confinement 
including crystals, plastic crystals, and liquid crystal (LC) phases.—"— LC phases in rod-like 
particles, for example, are observed in many natural and anthropogenic systems. Exam¬ 
ples include suspensions of colloidal particles such as vandium pentoxide (V 205 )^, and 
Gibbsite (A1(0H)3)^^, carbon nanotubes^^, and some biological systems such as protein 
fibers^^, tobacco mosaic virus^I, yd-virus^"— , polypeptide solutions^i^, and DNA— . Over 
the years, simple but non-trivial hard-core models have been used to study the formation 
of LC phases.— These models have played an important role to understand the behaviour 
of real systems. In particular, the hard-spherocylinder (HSC) has been used as a standard 
model to describe the LC behaviour of rod-like colloidal particles. The HSC model consist 
of a cylinder of length L and diameter D capped at each end by a hemisphere of diameter D, 
and it is shown in Fig. [I](a). Depending of the aspect ratio of the model, corresponding to 
the ratio L/D, suspensions of HSCs can exhibit the formation of isotropic, nematic, smec¬ 
tic, and solid phases. This rich phase behaviour has been confirmed by extensive computer 
simulat ions .— 

Despite our profound knowledge on the phase behaviour of rod-like particles, our un¬ 
derstanding of the phase behaviour of mixtures of anisotropic colloidal particles is still 
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limited due to the large parameter space that has to be explored, i.e., different combina¬ 
tions of concentrations, shapes, and sizes of the components, as well as different thermo¬ 
dynamic conditions. Experimental studies for mixtures of rod-like and spherical colloids 
have been reported.—^— From the modelling perspective, binary mixtures of HSC particles 
have been studied using computer simulations including rod-rod^*^*^!^, rod-diso^^— and 
rod-sphere^>^'— systems. These mixtures have shown the possibility of forming new struc¬ 
tures with properties which are difficult to attain in pure component mixtures. Rod-sphere 
mixtures are of particular interest as this case corresponds to one of the simplest colloid 
mixtures models, and the additional possibility of purely entropic depletion interactions, 
which can give rise to rich phase phenomena depending on the relative size ratio between 
the rods and the spheres.— 

The first statistical-mechanical theory to describe the isotropic-nematic phase transition 
of liquid crystal models was developed by Onsager. In his seminal work, Onsager— de¬ 
rived his simple density functional theory (DFT) for the isotropic-nematic transition by 
truncating the virial expansion at the level of second virial coefficient. The equilibrium 
state can then be determined by functional variation of the free energy with respect to the 
orientational distribution function. Although Onsager’s description is shown to be exact 
when the rods become infinitely long (because higher-order virial coefficients become neg¬ 
ligible decaying as D/Lr^), the theory does not accurately describe the phase behaviour of 
rod-like systems of intermediate values oi L/D when higher-order virial contributions are 
neglected. Several attempts have been made to extend Onsager’s theory by including the 
higher-order interactions. Recent progress in DFT— can provide appropriate approaches 
to the predictions of the thermodynamic properties of anisotropic fluids. A new free energy 
functional for inhomogeneous anisotropic fluids of arbitrary shape have been proposed within 
the framework of fundamental-measure theory^ which is based upon careful analysis of the 
geometry of the particles. Alternatively, the Parsons-Lee^“— approach provides a simple 
yet efficient way to incorporate the higher-order virial contributions which is neglected in 
Onsager’s method. Parsons^ proposed an approximation to decouple the orientational and 
translational degrees of freedom by mapping the properties of the rods to those of a refer¬ 
ence HS system. Lee^*^ approached the problem in a different way by introducing a scaling 
relation between virial coefficients of anisotropic particles and HS reference. Following two 
separate routes. Parsons and Lee reached the same expression for the free energy functional 
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which is commonly known as the Parsons-Lee (PL) theory2id2iiIii2ii2^. A straightforward 
extension of PL theory^ to the mixtures is the one-fluid approximation whereby one maps 
the mixture on to an effective one-component HS system. A decoupling approximation is 
used in the PL approach in which the system is represented as the effective hard sphere of 
the same diameter while any information about the geometry of the LC particles is included 
in the term of the factorized excluded volumes. In order to improve the PL treatment for 
mixtures, a many-fluid (MF) approach has been proposed^ where each component in the 
mixtures are mapped on to the corresponding effective HS system separately, thus LC mix¬ 
tures are represented as mixtures of HS. Following the separate routes of Parsons and of 
Lee, two versions of many-fluid theories can be developed: many-fluid Parsons (MFP) and 
many-fluid Lee (MFL) as alternatives for more accurate descriptions of LC mixtures. These 
many-fluid approaches have been assessed for a mixture of hard Gaussian particles and it 
has been shown that MFP is superior to the PL and MFL methods at moderate and high 
densities^. 

The focus of our current work is the isotropic-nematic phase behaviour of a HSC-HS 
mixture. Previous reports of the ordering in the HSC-HS binary system have been pre¬ 
sented including direct simulatioii^>^‘2^ and theoretical^i^SiiSS studies. The work of Cuetos 
and co-workers is of particular relevance: the one-fluid PL approach^ was used to study 
the isotropic-nematic phase diagram of the HSC-HS system characterized by rods of vari¬ 
ous lengths and diameters; comparisons where made with NPT Monte Carlo simulations^ 
employed to investigate the phase diagram and fluid structure of the mixtures. It is worth 
noting that in the NPT ensemble the system composition remains constant overall, which 
will lead to an inadequate description of the phase boundary as one enters metastable states 
which would otherwise phase separate into phases of distinct compositions. 

The purpose of our current work is twofold. First, the many-fluid Parsons theory is used 
to describe the HS-HSC system and comparisons are made with the one-fluid PL approach. 
It should be noted that in a one-component case both approaches reduce to the standard PL 
theory. Second, we present new Monte Carlo simulation results for the mixture. The local 
density (packing fraction), local composition, and orientational distributions are determined 
during the simulations to estimate the locations of the isotropic-nematic transitions of the 
mixture at various compositions in order to make a proper test of the accuracy of the two 
theories. 
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II. THEORY OF NEMATIC PHASE IN MIXTURES OF HARD 
PARTICLES 


In this section, the main steps leading to a formulation with both one-fluid and many- 
fluid theories are briefly recalled; further details can be found in Ref. [Sy. Consider an n- 


component mixture system of N unaxial (cylindrically symmetrical) hard anisotropic bodies 
in a volume U at a temperature T. The free energy functional of the system can be expressed 
as a contribution from an ideal (entropy) term (E‘^) and a residual (conhgurational) part 
(F’"®®): 


(3F _ 

'V ~ ~v~ ^ V 


( 1 ) 


where (3 = l/{k-QT) and k-Q is the Boltzmann constant; the temperature plays a trivial role in 
this case since only hard repulsive interactions between particles are considered. The ideal 
free energy accounts for the translational and orientational entropy and can be written as 


2=1 

where pi = Ni/V {N = number density of component i, p = N/V = 

J2i=iPiy Vi is the de Broglie volume of each species, incorporating the translational 
and rotational kinetic contributions of the ideal isotropic state. With the introduction of 
single-particle orientational distribution function fi{oj), the orientational entropy term cr[/j] 
can be expressed as an integration over all orientations cJ of a single particle: 


= J dufi{cj) ln[47r/i(cJ)]. 


( 3 ) 


For the residual part, Onsager’s original expressioi>^ is equivalent to truncating the virial 
expansion at second-virial level. At higher densities, however, the many-body correlations 
become progressively more and more important. Following the Parsons approach^, we can 
include higher-body contributions in an approximate manner. Assuming a pairwise additive 
hard interaction Uij{rki, 0 Jk,uji) between particle k of i-th component and particle I of j-th 
component, the pressure of a fluid mixture of n components can be written in the virial form 
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where I ^ k is used to avoid self interactions and (■ ■ •) represents the ensemble average. In 
the canonical ensemble, 


P = pk-BT - 



exp(-/317) (5) 


where the configurational partition function Z is dehned as 



( 6 ) 


and the total conhgurational energy is given by 



n n Ni Nj 


(7) 


i=i j=i k=i 

l^k 


The canonical pair distribution function is dehned as 


5'p(n2,cJ'i,CU2) 








da;^-2exp(-/3f/(r^%a;^^)), (8) 


where 5ij is the Kronecker delta. On integrating Equation (jS]), the expression for pressure 
can be written in a compact form: 






(9) 


Following the Parsons approach, the interparticle separation ri 2 is given in terms of the 
contact distance crij(ri 2 ,cJi,CJ 2 ) by dehning a scaled distance yij = 'ri 2 /(Jij(ri 2 ,cJi,cJ 2 ). The 
scaled distance does not explicitly depend on the orientations of the two particles and yij = 1 
corresponds to contact value. Using the dehnition of y^j, the pair distribution function (cf. 
Equation (jHj)) can be expressed as a function of scaled distance y^, i.e., pij = gij{y), which 
decouples the positional and orientational dependencies. In this way, a complicated pair 
potential Uij is mapped onto the spherically symmetrical hard-sphere potential: 




00 y < 1 


( 10 ) 
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and the expression for the pressure becomes 


1 /* ^ 
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( 11 ) 


where the excluded volume between a pair of particles is 0 J 2 ) = | / dfi 2 crfj(fi 2 , cJi, 0 / 2 ) 

and fi 2 = ri 2 /ri 2 . The form of hard repulsive pair interaction is a step function (cf. Equa¬ 
tion flTU]) h thus jdduij/dpij = — exp(/3Mjj)(5(|/jj — 1) (for example, see Ref. 1921) . Integrating 
over the scaled variable yij and noting that Ui+{y) = 0 when ?/ = 1, we then obtain 


^ n n f f 

P = pksT + - 

i=i j=i y y 


du2fi{(^i)fj{(^2)VP''{ufi, 0 / 2 ), ( 12 ) 


where gij{l^) ~ 5 fh®(l+) has been approximated as the corresponding hard-sphere contact 
value of pair distribution function. 

The residual free energy can then be obtained from the formal thermodynamic dehnition 
{dF/dV)NT = —P by integrating Equation flT^ over the volume: 


/3F’' 


V 


'^^'^^PiPjGij / dcJi / dU2fii(yJl)fj{(yj2)VP^{(Ji,(J2), 


(13) 


i=l j=l 

where Gij = p~^ dp'g^^(l~^). Onsager’s second-virial theory can be recovered with Gij = 1 

(i.e., 5 fh®(l’'') = 1) corresponding to the low-density limit. 

In this way, the theory of Parsons for a one-component fluid can be reformulated to 
describe a n-component mixtures of anisotropic bodies. As shown in Ref. [Sg, the standard 
“one-fluid” approach^ corresponds to Gij = G^^, where G^^ = p~^ dp'gQg(l~^), given in 
terms of the Carnahan-Starling form of the radial distribution function at contact^i^i^ 


. ( 14 ) 

with p = X)r=i and Vmj is the volume of the i-th species. The PL residual free energy 

can then be expressed as 
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Alternatively, in developing the many-fluid theory proposed in Ref. 
(volume) of each species individually, i.e.. 
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An expression for the contact value of the distribution function for hard-sphere mixture is 
given by Boublik— : 


HSjMix/i +\ 
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where the moments of the density are defined as (a = /Q) ^27=1 ^ ~ 0,1, 2, 3. 

Combining Equations flTdll and ffTTjl and noting the separate definition of Gij for each i — j 
pair, one obtains the many-fluid Parsons (MFP) form of the residual free energy irres,MFP_ 
In the one-component limit the contact value of the radial distribution function of 
the HS mixture (Equation (IT7l) ) reduces to the Carnahan-Starling expression (cf. Equa- 
tion flT^ ). Thereby, the MFP approach yields same descriptions as PL treatment for the 
pure-component systems. In the standard extension of the PL theory to mixtures one 
therefore adopts a van der Waals one-fluid (VDWl) approximation using an equivalent 
hard-sphere system with the effective diameter given by the VDWl mixing rule to represent 
the anisotropic mixtures. In contrast to the PL approach, each component is represented 
as a separate effective hard-sphere component, so that the excluded volume between a pair 
of i-th component and j-th component is weighted by the corresponding contact value of 
the HS mixture, . The equilibrium free energy of the system is determined from a 
functional variation with respect to the orientational distribution function fi{uj) of each 
component which leads to an integral equation for fi{u). The set of integral equations are 


solved numerically using an iterative procedure, details of which can be found in Ref. 
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In this work, we assess the adequacy of many-body theories such as the MFP for a binary 
mixture of hard spheres and hard spherocylinders. The models are depicted in Figure [T] the 
aspect ratio of the HSC is L/D = 5 and the diameter of the HS is taken to be the same as 
the diameter of the HSC, i.e., a = D. 

The excluded volumes corresponding to the HSC-HSC, HSC-HS and HS-HS interactions 
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FIG. 1. The hard-core models: (a) hard spherocylinder (HSC) of length L and diameter D; and 
(b) hard sphere (HS) of diameter a. In the current study, the length of the HSC is fixed to L = 5D 
and the diameter of the HS is the same value as that of the HSC, i.e., a = D. 

are given as 

Vi^rc-HSC = + 2 L^D \sinyl 

^HS-HS = (18) 

where 7 = arccos(a;i • U 2 ) is the relative orientation of the two HSC particles. The total 
excluded volume of the mixture is = x^gcHj|^c-Hsc + 2 xhsc 2 :hs 1 hsc-hs + ^hsIrs-hs 
where xhs and xhsc are mole fractions of HS and HSC species, respectively. Since the 
HSC particles are the anisotropic component in the system, /(cJ) is used to describe the 
orientation distribution of the HSC rods which is related to the nematic order parameter S 2 
of the system through 

82 = J duf{u) . (19) 

In particular, 5*2 = 0 corresponds to the isotropic state and S '2 = 1 for a perfectly-aligned 
nematic phase. 


III. MONTE CARLO SIMULATION OF PHASE COEXISTENCE IN 
MIXTURES OF HARD SPHERES AND HARD SPHEROCYLINDERS 

There are two common approaches to studying fluid-phase separation by molecular sim¬ 
ulation. Within the direct procedure the two coexisting phases are treated simultaneously 
in the presence of an interface with the usually periodic boundary conditions^^i^. The sta¬ 
bilization of a fluid interface corresponding to a system with a nonuniform density within a 
single simulation box is straightforward to implement with either molecular dynamics (MD) 
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or Monte Carlo (MC) techniques. This was first demonstrated by Croxton and Terrier— 
who performed MD simulations of the vapor-liquid interface of a Lennard-Jones system in 
two dimensions, and shortly afterwards by Leamy et ah— who stabilized the interface of 
a three-dimensional lattice gas (Ising model) by MC simulations. For a system which is 
sufficiently large (in the direction normal to the interface) one can simultaneously examine 
the bulk properties, in the central region of the coexisting phases as well as the interfacial 
properties. 

The direct molecular simulation of the isotropic-nematic (TN) phase transition in mix¬ 
tures of hard spheres and hard sp hero cylinders is particularly challenging because of the 
very low interfacial tension between the two phases; for example, the I-N interfacial ten¬ 
sion of a hard-core system of for thin disc-like particles has be estimated to be a few tenth 
of k-^T in units of the particle’s area^. As a consequence there is a very low energetic 
penalty associated with the deformation of the interface in such systems leading to large 
interfacial fluctuations; moreover, in the absence of an external held there is no resistance 
to the translation of a planar interface. The location of the bulk coexistence regions and 
the determination of the density and compositional prohles becomes a difficult task as a 
result. In order to break the symmetry of the system and reduce the effect of the interfacial 
huctuations one can introduce an external held by placing the system within structureless 
hard walls; this corresponds to removing the periodicity in one dimension (say the 2 ; direc¬ 
tion). An issue with this type of approach is that large systems have to be considered in 
order to study the true bulk phase behaviour and avoid capillary ehects. By keeping the 
separation between the walls large compared to the dimensions of the particles, one can 
simulate the phase coexistence in the hard-core HS-HSC mixtures with minimal ehect from 
the hard surfaces. 

Alternatively, the phase behaviour can be simulated using a popular Gibbs ensemble^^i^I 
in which the coexisting phases are retained in separate boxes and coupled volume changes 
and particles exchanges between the boxes are undertaken to meet the requirements of 
mechanical and chemical equilibria. However, in the case of hard anisometric particles, the 
acceptance ratio for the insertion of anisotropic particles will be very low, particularly at 
the high densities of the dense anisotropic phases of interest, requiring an impracticably 
large number of trial insertions for a proper equilibration of the systen>l^. A conventional 
simulation of the system within a single box will partially solve the problem since trial 
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insertions of the particles are no longer required. There is however a complication with the 
simulation of bulk phase equilibria of mixtures with a single simulation cell; though the 
phase transition between the various states can be traced as for a pure component system, 
the overall composition remains hxed preventing the fractionation of the different species in 
the various phases. 

In view of the aforementioned issues, we employ a less conventional NPnAT ensemble 
within a single cell where the component Pn of the pressure tensor normal to the interface 
is kept constant, so that the condition of mechanical equilibria is satished within the entire 
simulation cel l^^do^do^ ^ fjpg advantage of simulating the phase separation of mixtures by 
simultaneously considering the coexisting phases and the interface in a single cell is that 
this will allow for inhomogeneities in both the density and the composition of the system. 
By introducing an external held such as a hard surface one is able to examine both the bulk 
and interfacial regions of mixtures of hard core particles without constraining the density or 
composition of the individual bulk phases. 

We perform constant normal-pressure Monte Carlo simulation [NPnAT-MC) for a system 
of A^hsc = 1482 HSC particles of the aspect ratio L/D = 5 where the number of hard spheres 
is varied depending on mole fraction of the binary mixture a;HS,tot = 1 V'hs/(1V'hs + l^Hsc)- 
In this system, the intermolecular potential between any two particles is restricted to a 
pure repulsion. As shown in Figure [21 the simulation cell is a rectangular box of dimension 
Lx = Ly = 25D (corresponding to a hxed surface area A in the x-y plane oi A = 6250“^) 
and Lz varies according to the value set for P„. The parallel hard walls are positioned 
ai z = 0 and z = and standard periodic boundary conditions are applied in x and y 
directions. Since a hxed normal pressure is imposed along the ^ axis, the system volume in 
our NPnAT-MC simulation is allowed to huctuate by scaling the length of the 2 ; axis which 
moves the walls closer together or farther apart, while the system dimensions of the x and 
y axes and the x-y surface area are kept hxed. 

The NPnAT-MC simulation of the HSC-HS mixture is performed for 5 x 10® cycles 
to equilibrate the system and 5 to 8 x 10® cycles to obtain the average properties. Each 
MC cycles consist of = A^hs + -^hsc attempts to displace and rotate (in the case of a 
HSC particle) randomly chosen particles and one trial volume change corresponding to a 
contraction or extension in the .2 direction. The breaking of symmetry caused by the hard 
walls leads to inhomogeneous positional, orientational, and compositional distributions of 
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FIG. 2. The NPnAT-yiC simulation cell: hard walls are placed along the z axis and the cell is 
divided into three large bins: two surface regions close to the hard walls and a bulk region in the 
central part of the cell. A fixed normal pressure is imposed and the dimension of the system is 
allowed to fluctuate in the z direction. In this example a mixture system of hard spherocylinders 
(purple rods) and hard spheres (green spheres) is depicted. 

the system along the 2 ;-axis, so that the thermodynamic and structnral properties have to 
be determined locally. Smooth density and composition profiles are reqnired to identify the 
nniform region in the centre of the box which correspond to the bnlk phase. In order to 
evalnate the packing fraction rji{z), composition a;Hs (^)5 and order parameter S 2 {z) profiles, 
the simnlation box is divided into several bins of eqnal width 5z in the z direction; nbm = 200 
bins are used to calculate the packing fraction profile rji{z) = pi{z)V^^i,i = HS,HSC, where 
the number density profile of the component i is obtained from 

p,{z)='p^ * = HS,HSC (20) 

and the local composition is then obtained as a;Hs(^) = Phs(^)/(phs(^) + Phsc(^))- The 
packing fraction r^b and composition a:Hs,b of the bulk phase is then determined to be the 
values corresponding to the uniform regions of the packing fraction and composition profiles 
in the centre of the simulation cell. 

The nematic order parameter profile S 2 is obtained by determining the local nematic 
order parameter tensor Q(j) in each bin j: 

- ^ (ScJj ® cUj - I) 

i=l 



Q(j) = 


\ 2iVHsc 
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where iVHscj is the number of HSC particles in the j-th bin and I is the unit tensor. On 
diagonalising the tensor Q(j), three eigenvalues can be obtained and the largest eigenvalue 
dehnes the local nematic order parameter S 2 {z) of the j-th bin. Special care is required 
in calculating the order parameter profile because of finite-size effects. It has been shown 
by Eppenga and FrenkeU^i that the value of the nematic order parameter depends on the 
number of particles considered and the error in local order parameter is ^/W Nnsc/nun)- 
If we use ribin = 200 which is the same number of bins used for density profile, there 
are on average only 7(~ 1482/200) rods in each bin and the expected error in 5*2 (^) of 
~ 0.367 is large. Richter and Gruhnl^ have employed a methodology to correct for finite-size 
effects by introducing a function which bridges S 2 ^n^^c^oo{z) and *S'2 ,Vhsc(^) ^ly correlating 
the simulation data. A more direct way^ to reduce the error in the local nematic order 
parameter is to examine a larger system; for example, in the case of a system of 14820 rods 
(an order of magnitude larger than the system studied here) and 200 bins reduces the error 
in S 2 {z) to ~ 0.11. However, the simulation of a system of this size is very computationally 
expensive. As in our current work the focus is the bulk phase behaviour not the interfacial 
region, we divide the system into 3 large bins: 2 surface regions adjacent to the hard walls 
and the bulk region (cf. Figure [2]). With nbin = 3 system-size error in S 2 {z) decreases 
to ~ 0.04 where there are now an average of ~ 500 rods in each bin. The value of 5*2(^) 
corresponding to the central region is then taken to represent the nematic order parameter 
of the bulk phase S' 2 ,b- The reduced normal pressure is defined as P* = PnD/k-QT. 

IV. RESULTS AND DISCUSSION 
A. Pure hard spherocylinders 

Prior to demonstrating our results for mixtures of HS and HSC particles, it is instructive 
to begin by studying a system of pure HSC particles with aspect ratio oi L/D = 5. As is 
well known, the simple HSC model of a mesogen exhibits isotropic, nematic, smectic-A, and 
solid phases as the density of the system is increased^'^'^'^'^'^AS^. As an preliminary 
assessment we demonstrate that the bulk isotropic-nematic transition for the L/D = b HSC 
system contained between well separated parallel hard surfaces simulated using our NP^AT- 
MC methodology is essentially unaffected by the external field. The bulk phase behaviour 
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for the homogeneous system obtained using conventional constant pressure NPT-MC with 
full three dimensional periodic boundary conditions^ is compared with the corresponding 
data obtained using the constant normal-pressure NPnAT-MC methodology for the system 
between parallel hard surfaces in Figure Eland corresponding simulation data is reported in 
Table E 

The predictions of the isotropic and nematic branches with the MFP theory^ is also 
shown in Figure [3] which reduces to the commonly employed PL theory^"— for one- 
component system: the theory is seen to provide a reasonably quantitative description 
of the isotropic and nematic branches of the equation of state and the position of the 
isotropic-nematic transition. We can infer that the results obtained for the HSC particles 
contained between the parallel hard surfaces are fully consistent with those of the fully 
periodic homogeneous system, confirming that at least for this system size and geometry of 
the simulation cell the presence of hard surfaces only has a small (stabilizing) effect on the 
isotropic-nematic transition; in this case the separation between the surfaces, which defines 
the dimension the box in the .2 direction, ranges from ~ 899.85D at the lowest density 
{P* = 0.001) studied to ~ 23.57/1 for the highest density {P* = 1.26) nematic state. 
Example of the packing fraction (density) profiles ri{z) obtained for a low-density isotropic 
state, an intermediate-density nematic state, and a moderately high-density nematic state 
are displayed in Figure HI the flat part of the prohles correspond to the homogenous bulk 
phases in the central region of the simulation cell which allows one to determine the equi¬ 
librium bulk density r]i, with confidence. 

As in other studies of confined liquid-crystalline systemsl2£Ai^^ signihcant structure is 
also apparent close to the surfaces; a full analysis of the surface effects such as wetting, 
de-wetting, surface nematization, and adsorption will be left for future work. 

In order to estimate the location of the bulk isotropic-nematic bulk phase transition, we 
examine the density dependence of the nematic order parameter in Figure |5l The order pa¬ 
rameter of a hnite system S2 ,Niisc (^) converges quickly to the limiting bulk thermodynamic 
value 5'2,AfHsc^oo(^) states with intermediate to high orientational order {S 2 ,Nu 5 c ~ O-^)- 
In Figure Owe display the order parameter profiles for states of low to moderate orienta¬ 
tional order (0.1 < 5'2,Whsc ~ 0 -^) close vicinity of the isotopic-nematic transition. 

Snapshots of typical configurations of these two states are shown in Figured the low-density 
conhguration corresponds to an bulk isotropic state, the intermediate-density configuration 
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TABLE I. Constant normal-pressure MC {NPnAT-MC) simulation results for bulk isotropic- 

nematic phase behaviour of Nhsc = 1482 hard spherocylinders with an aspect ratio L/D = 5. 

The reduced normal pressure P* is set in the simulation and corresponding bulk values of packing 

fraction nematic order parameter S' 2 ,b) and box length are obtained as configurational av¬ 

erages. The isotropic phase is denoted by Iso, the nematic by Nem, and the pre-transitional states 

by Pre. 

P* 77b *S' 2 ,b L^jD Phase 

0.001 

0.004 


0.042 

899.85 

Iso 

0.003 

0.011 


0.042 

881.84 

Iso 

0.005 

0.018 


0.042 

562.67 

Iso 

0.01 

0.031 


0.042 

327.44 

Iso 

0.02 

0.057 


0.043 

183.67 

Iso 

0.05 

0.098 


0.045 

105.10 

Iso 

0.10 

0.144 


0.046 

71.12 

Iso 

0.20 

0.199 


0.049 

51.17 

Iso 

0.30 

0.235 


0.053 

43.39 

Iso 

0.40 

0.271 


0.056 

38.84 

Iso 

0.50 

0.300 


0.067 

35.12 

Iso 

0.60 

0.318 


0.072 

32.81 

Iso 

0.70 

0.340 


0.074 

30.85 

Iso 

0.80 

0.356 


0.085 

29.23 

Iso 

0.90 

0.372 


0.099 

27.68 

Iso 

1.00 

0.386 


0.116 

26.75 

Iso 

1.10 

0.397 


0.371 

25.40 

Pre 

1.12 

0.401 


0.456 

25.18 

Pre 

1.14 

0.419 


0.553 

24.88 

Nem 

1.16 

0.420 


0.555 

24.82 

Nem 

1.18 

0.425 


0.570 

24.61 

Nem 

1.20 

0.433 


0.619 

24.22 

Nem 

1.22 

0.434 


0.692 

23.95 

Nem 

1.24 

0.435 


0.727 

23.73 

Nem 

1.26 

0.443 


0.f^9 

23.57 

Nem 






Vb 


FIG. 3. The isotropic-nematic phase behaviour of hard-spherocylinder (HSC) rods with an as¬ 
pect ratio of L/D = 5. The simulation data obtained for the system of Nhsc = 1482 particles 
contained between parallel hard surfaces using our NPnAT-M.C approach (filled symbols), where 
P* = PnD^ / (k^T) is the dimensionless normal pressure and r]y^ = Pbl^nsc is the bulk value of the 
packing fraction in the central region of the cell, are compared with the corresponding simulation 
data obtained by McGrother et ai— for the fully periodic system using the conventional NPT-M.C 
approach (open symbols), where now P* = P/{kB,T) and r] = pVnsc are the values for the homo¬ 
geneous system. The circles corresponds to the isotropic states, the triangles to the nematic state, 
and the squares to the pre-transitions states. The curve is the theoretical predictions obtained with 
the PL theory. An enlargement of the isotropic-nematic transition region is shown in the inset. 

to a pre-transitional state, while the high-density conhguration has clearly undergone a tran¬ 
sition to a bulk nematic phase. The pre-transition states, assumed here to correspond to 
nematic order parameters in the range 0.3 < 5'2,AfHsc ~ system size effects 

and can also be exacerbated by the conhnement and potentially very slow relaxation order 
processes in the form of slow nucleation kinetics. 

The discrepancy between the values of the nematic order parameter obtained for prohles 
with Ubin = 3 and ribin = 20 histogram bins in the z direction is very small so that the 
size effects are expected to be small in this case. An isotropic bulk phase region is clearly 
found in the case of the state with a pressure of P* = 1.00; in this case the orientational 
order in the bulk region is found to be low corresponding to an nematic order parameter 
of 'S' 2,6 = 0.116. The order parameter prohle for the pre-transitional state with a normal 
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FIG. 4. Packing fraction profile rj{z) for Nhsc = 1482 pure hard-spherocylinder (HSC) rods with 
normal pressure P* = 0.2 to 1.2 obtained using NPnAT-M.C. Two hard walls are positioned at 
z = 0 and z = L^, where is the length of the z axis. 

pressure of = 1.12 exhibits a curve with no uniform bulk region; the average of the 
nematic order parameter of 5'2,b = 0.456 obtained in the central part of the cell does not 
therefore represent that of a true bulk phase. In the case of the denser state corresponding 
to P* = 1.14, the value of the bulk nematic order parameter 5*2,b = 0.553 obtained as an 
average over Ubm = 3 bins is essentially equivalent to the value of 5*2,b = 0.556 obtained 
with TT-bin = 20 bins in the homogeneous central region of the simulation cell. The difference 
in the orientational ordering for the states corresponding to normal pressures of P,^ = 1.12 
and P* = 1.14 is also apparent from Figure [3 where the orientations of HSC rods have 
been colour coded to aid the visualization. Small clusters of nematic domains are seen for 
the pre-transitional states (Fig[6l), which lead to slightly larger values of the nematic order 
parameter S 2 ,b 

Clearly, the equilibrium state at the pressure of P^ = 1.14 corresponding to a bulk density 
of r]f, = 0.419 and nematic order parameter of 5*2,b = 0.553 can be taken to correspond to the 
lowest-density nematic state, while the state at the slightly lower pressure of P^ = 1.12 is seen 
to exhibit some small nematic clusters which would correspond to a pre-transitional state; 
large region with random orientations corresponding to a bulk isotropic liquid can be seen 
in the case of the system with P* = 1.00. Our estimate of the isotropic-nematic transition 
for the L/D = b HSC system from an analysis of this data is ?7b,iso = 0.386, r7b,nem = 0.419 
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FIG. 5. The density dependence of the bulk nematic order parameter S' 2 ^b for hard-spherocylinder 
(HSC) rods with aspect ratio L/D = 5. The simulation data obtained for the system of iV = 1482 
particles contained between parallel hard surfaces using our NPnAT-M.C approach (filled symbols), 
where 5'2,b is the bulk value of the nematic order parameter and r/b = /Obk^HSC is the bulk value of the 
packing fraction in the central region of the cell, are compared with the corresponding simulation 
data obtained by McGrother et ah^ for the fully periodic system using the conventional NPT- 
MC approach (open symbols), where now S2 and r] = pVusc a-re the values for the homogeneous 
system. The circles correspond to the isotropic states, the triangles to the nematic states, and the 
squares the pre-transitional states in the vicinity of isotropic-nematic transition. The curve is the 
theoretical prediction obtained with the PL theory. 

for the bulk coexisting phases which are in good agreement with the corresponding re- 
snlts estimated from conventional 7VPT-MC simnlation with fnll three-dimensional periodic 
bonndary conditions (r^iso = 0.407, ?7nem = 0.415)^; the coexistence pressnre is arithmetic 
average between the normal pressures corresponding to the highest-density isotropic state 
and the lowest-density nematic state, P* = 1.07, which is slightly lower than that estimated 
for the system with full periodicity {P* = 1.19)^. A determination of the free energy and 
chemical potential of the system would allow one to get a more precise estimate of the po¬ 
sition of the phase transition, but this is beyond the scope of the current study. A slight 
stabilization of the isotropic-nematic transition (corresponding to a lowering the transition 
packing fraction by 1 to 2%) is therefore found in our systems of HSC rods placed between 
two parallel hard walls. One would certainly expect surface induced nematization to sta- 
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FIG. 6. The nematic order parameter S 2 {z) (bottom panel) and packing fraction q(z) (top panel) 
profiles for hard-spherocylinder (HSC) rods with an aspect ratio oi L/D = 5 obtained by simulating 
the system of = 1482 particles between parallel hard surfaces placed along the z direction. 
Isotropic {P* = 1.00), pre-transitional {P* = 1.12), and nematic {P* = 1.14) states in the vicinity 
of the isotropic-nematic transition are examined; the thermodynamic state is characterized by the 
value of the normal pressure tensor, P* = P^D^/{k^T). In the case of S 2 {z) the profiles are 
constructed from histograms with nbin = 3 (symbols) and Ubin = 20 (curves) to assess possible 
system size effects. 

bilize the transition to a bulk nematic phase. These surface effects are however not the 
focus of the current study and will be discussed in detail in subsequent studies. When the 
two conhning walls are well separated the effect on the bulk isotropic-nematic transition is 
inappreciable. However, the effects of conhnement will become increasingly more important 
as the surfaces are brought closer together, i.e., at higher pressures and/or for small system 
sizes. For the systems studied in our current work the two hard walls in our simulation cells 
are far enough apart (even for the highest density states) so that the effect of the surfaces 
on the bulk phase transition is small. 


B. HSC-HS mixture 

We now turn our attention to a binary mixture of HS and HSC particles characterized by 
overall HS composition xhs ranging from 0.05 to 0.20. As with the pure component system 
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FIG. 7. Snapshots of typical configurations of Nhsc = 1482 hard-spherocylinder (HSC) rods 
with an aspect ratio of L/D = 5 obtained by simulating the system between parallel hard surfaces 
places along the z axis, (a) Isotropic {P* = 1.00), (b) pre-transitional {P* = 1.12), and (c) nematic 
{P* = 1.14) states in the vicinity of the isotropic-nematic transition are examined (see the Caption 
of Figure[6]for further details). The colours quantify the orientation of the rod-like particles relative 
to the frame of the simulation box. 


the mixture is contained between parallel hard structureless surfaces and the isotropic- 
nematic phase transition is simulated using the NP^AT-ViC: method (c/. Section HTT)) . As 
well as being of different density the coexisting isotropic and nematic phases will exhibit a 
fractionation of the two components, with an accumulation of the the rod-like particles in 
the orientationally ordered nematic phase. The fluid phase separation in hard-core system 
of this type is an entropy driven process. This is not to be confused with the depletion 
driven phase behaviour exhibited by anisometric colloidal particles on addition of polymer 
where the polymer induces an effective attractive interaction (depletion force) between the 
colloids, that would otherwise only interact in a purely repulsive fashion, 
van der Waals like “vapour-liquid” transition (see for example references 


58 

66 

68 

69 

114 


llTj and the excellent monograph by Lekkerkerker and Tuinier— ). The simulated phase 


boundaries of our HSC-HS mixture will be compared with the theoretical corresponding 
predictions obtained with the one-fluid PL and two-fluid MFP approaches. A simulation 
cell in a low-density thermodynamic state is slowly compressed and equilibrated to obtain the 
dependence of bulk packing fraction and bulk composition on the equilibrium bulk pressure 
(which for our system with planar symmetry also corresponds to the normal component of 
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the pressure tensor, P„). The bulk values of the phases are again obtained as averages of 
the density and composition profiles in the central region of the simulation cell. 

Typical density, composition and nematic order parameter profiles for bulk isotropic, 
intermediate isotropic-nematic, and bulk nematic states of mixtures with overall HS compo¬ 
sition xhs = 0.10 are shown in Figure [HI respectively. As in the case of the pure-component 
HSC system, the density prohles of the HSC-HS mixture reveal significant order of the HSC 
close to the walls and a well-defined bulk region characterised by the flat profiles in the 
central portion of the simulation cell. This positional and orientational order is further con- 
hrmed by the large values of composition and nematic order parameter of the HSCs, and 
the low concentration of HSs near the surface. The characteristic flat profiles of the nematic 
order parameter in the bulk isotropic at P* = 1.10 and nematic at P* = 1.13 states can be 
seen in Fig. [H](c); the V-shape nematic order profile of the pre-transitional state at P* = 1.21 
is also clearly apparent. The corresponding data for the mixture with 10% HS is given in 
Table in The dependence of the bulk packing fraction r^b on the applied normal pressure P* 
for the HSC-HS system with an overall composition of xns.tot = 0.1 is illustrated in Figure 
[9]^a). The theoretical description of isotropic and nematic branches of the equation of state 
obtained with the PL and MFP approach at the same overall composition are also plotted in 
Figure [Oja) for comparison. There is only a negligible difference between the PL and MFP 
results in the isotropic state. This is to be expected, since both the density and the ordering 
in this region is small and as a consequence of the effective hard-sphere treatment of the 
excluded volume contributions with a one- or two-fluid approximation should be similar for 
the isotropic state. One should note that while within the theories (both PL and MFP) the 
density and composition are input variables for a given system and the pressure is an out¬ 
put, while for our NPnAT-MC simulation the equilibrium pressure is specihed and the bulk 
density and composition are obtained as averages from the central region of the cell. An ob¬ 
servable difference between the one- and two-fluid theories can be seen in the vicinity of the 
isotropic-nematic transition where the branches of the equation of state obtained by simula¬ 
tion data experience an abrupt change in slope indicating the transition to the orientationally 
ordered state. From the results depicted in Figure [HI one can infer that predictions with the 
two-fluid MFP approach is marginally superior to that with the one-fluid PL at least for 
the nematic phase. The isotropic-nematic transition point can be identihed from the abrupt 
change in nematic order parameter as is apparent from Figure [9jb). The typical snapshots 
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of configurations for the highest-density bulk isotropic state (r^b = 0.414, XHs,b = 0.109) and 
the lowest-density bulk nematic state {rjh = 0.419, a;HS,b = 0.108) are also included in order 
to visualise differing degrees of orientational order. 

TABLE II. Constant normal-pressure MC (NPnAT-MC) simulation results for bulk isotropic- 
nematic phase behaviour of mixtures of A^hsc = 1482 HSC and A^hs = 165 HS particles for an 
overall HS composition of xns.tot = 0.10. The reduced normal pressure P* is set in the simulation 
and corresponding bulk values of packing fraction r]y^, composition a;HS,b) nematic order parameter 
S' 2 ,b) and box length L^- are obtained as configurational averages. The isotropic phase is denoted 
by Iso, the nematic by Nem, and the pre-transitional states by Pre. 


p* 

hb 

a^HS.b 

52,b 

L,/T» 

Phase 

0.70 

0.336 

0.104 

0.069 

32.07 

iso 

0.80 

0.355 

0.105 

0.076 

30.34 

iso 

0.90 

0.370 

0.107 

0.092 

28.58 

iso 

1.00 

0.383 

0.109 

0.118 

27.56 

iso 

1.10 

0.396 

0.109 

0.167 

26.71 

iso 

1.14 

0.411 

0.109 

0.239 

26.15 

iso 

1.16 

0.414 

0.109 

0.286 

25.84 

iso 

1.18 

0.408 

0.115 

0.317 

25.65 

pre 

1.20 

0.413 

0.114 

0.399 

25.3 

pre 

1.21 

0.414 

0.114 

0.345 

25.35 

pre 

1.22 

0.416 

0.114 

0.339 

25.27 

pre 

1.23 

0.419 

0.108 

0.551 

25.18 

nem 

1.24 

0.420 

0.108 

0.540 

24.97 

nem 

1.25 

0.421 

0.112 

0.584 

24.66 

nem 

1.26 

0.423 

0.121 

0.614 

24.58 

nem 

1.28 

0.430 

0.111 

0.661 

24.35 

nem 

1.30 

0.436 

0.11 

0.708 

24.13 

nem 

1.32 

0.441 

0.109 

0.724 

23.89 

nem 


The corresponding results for the HSC-HS system with the highest overall HS composition 
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FIG. 8. (a) The packing fraction rji{z), (b) composition x\{z), and (c) nematic order parameter 
5*2( 2 ;) profiles for mixtures of A^hsc = 1482 hard-spherocylinder (HSC) and A^hs = 165 hard-sphere 
(HS) particles for an overall composition of xns.tot = 0.10. Typical bulk isotropic {P* = 1.10), 
pre-transitional {P* = 1.21), and nematic {P* = 1.23) states are examined. The packing fraction 
and composition profiles profiles of HS particl^^are shown in the insets of (a) and (b). 








































TABLE III. Constant normal-pressure MC {NPnAT-M.C) simulation results for bulk isotropic- 
nematic phase behaviour of mixtures of Ahsc = 1482 HSC and A^hs = 371 HS particles for an 
overall HS composition of xns.tot = 0.20. The reduced normal pressure P* is set in the simulation 
and corresponding bulk values of packing fraction composition xns.b; nematic order parameter 
iS' 2 ,b) and box length are obtained as configurational averages. The isotropic phase is denoted 
by Iso, the nematic by Nem, and the pre-transitional states by Pre. 


p* 

hb 

aJHS.b 

52,b 

LJD 

Phase 

0.80 

0.348 

0.206 

0.131 

31.05 

Iso 

0.90 

0.361 

0.211 

0.151 

29.75 

Iso 

1.00 

0.379 

0.212 

0.172 

28.68 

Iso 

1.10 

0.391 

0.213 

0.163 

27.85 

Iso 

1.20 

0.403 

0.222 

0.192 

26.45 

Iso 

1.22 

0.403 

0.224 

0.193 

26.21 

Iso 

1.24 

0.404 

0.220 

0.260 

25.95 

Iso 

1.26 

0.408 

0.226 

0.301 

25.62 

Pre 

1.28 

0.414 

0.225 

0.404 

25.53 

Pre 

1.30 

0.415 

0.230 

0.415 

25.35 

Pre 

1.31 

0.417 

0.233 

0.419 

25.26 

Pre 

1.32 

0.430 

0.219 

0.619 

25.02 

Nem 

1.34 

0.430 

0.223 

0.612 

24.71 

Nem 

1.36 

0.431 

0.226 

0.628 

24.50 

Nem 

1.38 

0.432 

0.226 

0.619 

24.64 

Nem 

1.39 

0.436 

0.226 

0.650 

24.53 

Nem 


2^HS,tot = 0.20, are shown in Figure [TO] and in Table HTTl The isotropic-nematic transition lies 
somewhere between highest-density bulk isotropic state at P* = 1.24 and the lowest-density 
bulk nematic state at P* = 1.32 where a clear change in the density and nematic order 
parameter is exhibited: the values of the packing fraction and nematic order parameter for 
these two states are r^b = 0.404 and S' 2,6 = 0.260 for the bulk isotropic phase and r^b = 0.430 
and S' 2,6 = 0.619 for the bulk nematic phase. For completeness the less extensive data for 
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overall hard-sphere compositions of 5% and 15% are given in Tables [I^ and |Vl 

The isotropic-nematic phase boundaries estimated for the HSC-HS mixtures for bulk 
phase compositions ranging from the pure HSC system {xnsfi = 0) to XHS,b ~ 0.3 is shown 
in Figure [11] Here, we also compare the predictions of PL and MFP theories with our 
NPnAT-MC simulations as well as with the previously reported iVPT-MC data for sys¬ 
tems in full three-dimensional periodic boundary conditions^. The one-fluid and two-fluid 
approaches are both seen to describe the coexisting packing fractions as monotonically in¬ 
creasing functions of the bulk composition for the isotropic and nematic phases. For the 
phase boundary of the nematic states found at the higher packing fractions, the difference 
between the description obtained with the PL and MFP approaches becomes more marked 
with increasing composition of the spherical particles. Our simulated values for the nematic 
and isotropic phase boundaries are consistent with those obtained with a fully periodic sys¬ 
tem; the simulation data are in reasonably good agreement with the PL and MFP theoretical 
descriptions. There is an overestimate of the hrst-order character of the isotropic-nematic 
transition with scaled Onsager theories of this type based on a underlying description at the 
level of the second virial coefficient^. The approximations inherent in mapping HSC-HS 
mixture on to an equivalent HS mixture in order to simplify the computation of higher- 
order virial contribution may also lead to an exaggeration of the hrst-order nature of the 
transition. On increasing the overall composition of the HS particles,the density gap be¬ 
tween the isotropic and nematic coexisting states is seen to become larger as found with 
our simulations. By contrast,the density gap between the phase boundaries obtained from 
conventional NPT-MC simulations of the fully periodic systent^ appears to shrink slightly. 
In fully periodic NPT-MC simulations of this type the system remains essentially homoge¬ 
neous so that the composition of the state remains hxed. As the pressure is increased the 
system will undergo a transition from an isotropic to nematic phase but the states will be 
constrained to have the same composition. As a consequence of lack of fractionation of the 
species between the two phases with the NPT-MC simulations it is difficult to differentiate 
metastable states within the binodal region from those corresponding to the coexistence 
boundaries. 

The NPnAT-MC simulation data for the HSC-HS mixture are also summarized reported 
in Table I VI I where the slight composition asymmetry between coexisting isotropic and ne¬ 
matic phases is clearly apparent. For the HSC-HS system, the addition of spherical particles 
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is found to destabilize. The formation of a bulk nematic state predominantly composed of 
HSC rods will cause a reduction in the concentration of the HS particles in the same region, 
and as a consequence the orientational ordered state which is of higher density but lower 
bulk HS composition coexists with an isotropic state of lower-density and higher bulk HS 
composition. 

Finally, a phase diagram in the pressure-composition {P*-xi{s,h) plane is shown in Figure 
[121 A very narrow region of isotropic-nematic coexistence is obtained with our NPnAT- 
MC simulation approach. The coexistence region is seen to be at lower pressures that 
that predicted with the PL and MFP theories or obtained by fully periodic iVPT-MC 
simulatioii^. It should be noted that the results obtained at higher compositions with the 
fully periodic simulations is in good agreement with the nematic branch predicted with MFP 
theory. The quality of both PL and MFP is affected by a shift in pressure for the pure HSC 
system. The predictions with the two-fluid MFP theory are seen to be much better than 
one-fluid PL theory, particularly for systems of higher composition. 


V. CONCLUSIONS 

In this paper we have studied phase behaviour of mixtures of purely repulsive rod-like 
(HSC) and spherical (HS) particles. Using the new NPnAT Monte Carlo simulations we have 
constructed the isotropic-nematic phase diagrams of the HSC-HS mixture in the pressure- 
density, pressure-concentration and density-concentration projections. The comparison of 
our results with previously reported fully periodic iVPT-MC data^ reveals good agreement 
at low HS concentrations but some discrepancy occurs for higher HS concentrations espe¬ 
cially in the nematic phase. In conventional iVPT-MC simulations the system is essentially 
homogeneous such that the overall composition of the system remains fixed. This means 
that unless one has a very large system the coexisting isotropic and nematic states would 
have identical composition a^^^d^o . By contrast, using our NPnAT-M.C algorithm we have ob¬ 
served compositional asymmetry between the coexisting phases, with a slight increase of the 
HSC concentration in the nematic phase due to packing effects. This conclusion is also sup¬ 
ported by the fact that our NPnAT results for pure HSC particles are completely consistent 
with the available data obtained from iVPT-MC simulations using full three-dimensional 
boundary conditions^. 
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FIG. 9. (a) The pressure-density P^-ilh phase diagram, and (b) nematic orientational order pa¬ 
rameter S 2 ^b of mixtures of hard spherocylinders (HSC) and hard spheres (HS). The results of 
NPnAT-MC simulations for mixtures of A^hsc = 1482 HSC and A^hs = 165 HS particles for an 
overall HS composition of XHS,tot = 0.10 contained between well separated parallel hard walls 
are represented as: circles for the isotropic states; triangles for the nematic states; and squares 
for the pre-transitional states. The continuous and dashed curves represent the predictions using 
the one-fluid PL and two-fluid MFP theories, respectively; the low-density branch corresponds to 
isotropic states and the high-density branch to nematic states. The snapshots in (b) correspond 
to the highest-density bulk isotropic state (77b = 0.414,3:Hs,b = 0.109) and the lowest-density bulk 
nematic state (% = 0.419, XHS,b = 0.108). 
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FIG. 10. (a) The pressure-density P^-Vh phase diagram and (b) nematic orientational order param¬ 
eter S 2 ,b of mixtures of hard spherocylinders (HSC) and hard spheres (HS) particles. The results 
of NPnAT-MC simulations for mixtures of A^hsc = 1482 HSC and A^hs = 371 HS particles for 
an overall HS composition of a:HS,tot = 0.20 contained between well separated parallel hard walls 
are represented as: circles for the isotropic states; triangles for the nematic states; and squares 
for the pre-transitional states. The continuous and dashed curves represent the predictions using 
the one-fluid PL and two-fluid MFP theories, respectively; the low-density branch corresponds to 
isotropic states and the high-density branch to nematic states. The snapshots in (b) correspond 
to the highest-density bulk isotropic state (r^b = 0.404, XHS,b = 0.220) and the lowest-density bulk 
nematic state (r^b = 0.430, XHs,b = 0.219). 
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FIG. 11. The isotropic-nematic phase diagram of mixtures of hard spherocylinders (HSC) and hard 
spheres (HS) particles in ?7b(?7)-XHs,b plane. The results of NPnAT-MC simulations for mixtures 
of A^hsc = 1482 HSC and A^hs = 165 HS particles for an overall HS composition of xns.tot = 0.1 
are represented as filled circles for the isotropic branch and filled triangles for the nematic branch. 
The NPnAT-MC data are compared with predictions obtained with the one-fluid PL (dashed 
curves) and two-fluid MFP (continuous curves) theories, and with previously reported NPT-MC 
simulation data for fully periodic systems: open squares^ and crosses^. 

The advantage of our method is that the fluid phase separation can be studied within 
a single simulation cell which circumvents the problem of inserting anisotropic particles 
inherent in other techniques such as Gibbs ensemble MC (GEMC). Additionally, particle 
exchanges are allowed between the surface and bulk regions so that one is able to treat the 
two bulk states coexisting at different bulk densities (packing fractions) as well as different 
bulk compositions. The effect of the pair of the auxiliary hard walls put at the end of the 
box is not appreciable in the bulk region for sufficiently prolonged systems in the direction 
perpendicular to the walls and considerably helps to stabilise the phase separated states 
with a low interfacial tension which as typically exhibited by hard-body systems. 

We further compared the new simulation data with two theoretical predictions that go 
beyond the Onsager second-virial theory and extend the well known PL theory to mixtures 
in two different ways: using the one-fluid approximation whereby one maps the mixture 
onto an effective HS mixture, and by treating both species separately as an effective HS 
mixture which we refer to as many-fluid approximation. The comparison reveals that both 
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^HS,b 

FIG. 12. The isotropic-nematic phase diagram of mixtures of hard spherocylinders (HSC) and hard 
spheres (HS) particles in plane. The results of NPnAT-MC simulations for mixtures of 

-^HSC = 1482 HSC and A^hs = 165 HS particles for an overall HS composition of xns.tot = 0.1 
are represented as filled circles for the isotropic branch and filled triangles for the nematic branch. 
The NPnAT-yiC data are compared with predictions obtained with the one-fluid PL (dashed 
curves) and two-fluid MFP (continuous curves) theories, and with previously reported NPT-M.C 
simulation data for fully periodic systems; open squares^ and crosses^. The simulation results 
obtained from NP^AT simulations have been rescaled with the pure component pressure in the 
inset; the same rescaling procedure applied to the results obtained using the PL and MFP theories. 

the one- and many-fluid approaches provide a reasonably accurate quantitative description 
of the mixture including the isotropic-nematic phase boundary and degree of orientational 
order of the HSC-HS mixtures. The many-fluid prediction of the coexisting pressure is 
arguably found to be slightly better to that obtained with the one-fluid method. However, 
systems with larger aspect ratios should be considered to make a better assessment of the 
performance of the theories. 

Our work can be directly extended in a number of ways. For instance, we have restricted 
our attention to systems of HSC rods with an aspect ratio of L/T) = 5 and HS particles 
of diameter a = D. Considering larger aspect ratios would be a more stringent test of the 
theoretical predictions, as one would expect a more signihcant improvement of the many- 
fluid treatment compared to the one-fluid approach. The compositional asymmetry in the 
coexisting isotropic and nematic phases is expected to be more appreciable. The method for 
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TABLE IV. Constant normal-pressure MC (NP^AT-MC) simulation results for bulk isotropic- 
nematic phase behaviour of mixtures of Vhsc = 1482 HSC and Vhs = 78 HS particles for an 
overall HS composition of xns.tot = 0.05. The reduced normal pressure P* is set in the simulation 
and corresponding bulk values of packing fraction r/b, composition XHS,b) nematic order parameter 
5 *2,b) and box length are obtained as configurational averages. The isotropic phase is denoted 
by Iso, the nematic by Nem, and the pre-transitional states by Pre. 


p* 

hb 

a^HS,b 

52,b 

L,/D 

Phase 

0.700 

0.337 

0.052 

0.072 

31.380 

Iso 

0.800 

0.355 

0.054 

0.080 

29.610 

Iso 

0.900 

0.372 

0.054 

0.100 

28.210 

Iso 

1.000 

0.382 

0.056 

0.123 

27.190 

Iso 

1.100 

0.401 

0.057 

0.245 

25.860 

Iso 

1.120 

0.405 

0.057 

0.298 

25.630 

Iso 

1.140 

0.409 

0.059 

0.345 

25.330 

Pre 

1.150 

0.416 

0.058 

0.392 

25.330 

Pre 

1.160 

0.412 

0.061 

0.392 

25.270 

Pre 

1.170 

0.418 

0.057 

0.574 

24.960 

Nem 

1.180 

0.424 

0.056 

0.657 

24.540 

Nem 

1.200 

0.429 

0.062 

0.641 

24.520 

Nem 


locating the phase boundaries can be improved using thermodynamic integration. Further¬ 
more, our NPnAT-MC approach can be directly applied to describe the surface phenomena 
due to both fluid-fluid and fluid-wall interfaces. These interfaces plays an important role 
in liquid crystalline system o^^^d^^ due of rich surface-induced effects (e.g., nematization and 
smectization), characteristic in systems comprising anisotropic particles. 
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TABLE V. Constant normal-pressure MC (NPnAT-MC) simulation results for bulk isotropic- 
nematic phase behaviour of mixtures of A^hsc = 1482 HSC and A^hs = 262 HS particles for an 
overall HS composition of xns.tot = 0.15. The reduced normal pressure P* is set in the simulation 
and corresponding bulk values of packing fraction r/b, composition XHS,b) nematic order parameter 
5 *2,b) and box length are obtained as configurational averages. The isotropic phase is denoted 
by Iso, the nematic by Nem, and the pre-transitional states by Pre. 


p* 

7b 

a;HS,b 

52,b 

L,ID 

Phase 

1.000 

0.381 

0.163 

0.151 

27.948 

Iso 

1.100 

0.392 

0.166 

0.180 

26.819 

Iso 

1.200 

0.412 

0.169 

0.260 

25.524 

Iso 

1.210 

0.406 

0.174 

0.241 

25.645 

Iso 

1.220 

0.412 

0.172 

0.306 

25.757 

Iso 

1.250 

0.421 

0.171 

0.360 

25.290 

Pre 

1.260 

0.420 

0.171 

0.441 

25.104 

Pre 

1.270 

0.428 

0.167 

0.612 

24.805 

Nem 

1.280 

0.424 

0.170 

0.560 

24.990 

Nem 
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TABLE VI. The isotropic-nematic transition estimated from NP^AT-MC simulations of mixtures 
of Vjjsc = 1482 hard spherocylinder (HSC) and Vhs hard-sphere (HS) particles for varying over¬ 
all HS compositions of XHS,tot contained between well separated parallel hard walls. The normal 
pressure P* = PnD^/(k-QT) is set during the simulation, bulk packing fractions rjh and bulk com¬ 
positions XHS,b of coexisting isotropic (iso) and nematic (nem) states are obtained as averages from 
the central region of the cell. The bulk nematic order parameter £' 2 ,b of the lowest-density nematic 
bulk phase is also shown. 


a^HS,tot 

p* 

rx,iso 

^b,iso 

^HS,b,iso 

p* 

'^n,nem 

^b,nem 

^HS,b,nem 

52,b 

0 

1.00 

0.386 

0 

1.14 

0.419 

0 

0.553 

0.05 

1.12 

0.405 

0.057 

1.17 

0.418 

0.057 

0.574 

0.10 

1.16 

0.414 

0.109 

1.23 

0.419 

0.108 

0.551 

0.15 

1.22 

0.412 

0.172 

1.27 

0.428 

0.167 

0.612 

0.20 

1.26 

0.408 

0.226 

1.32 

0.430 

0.219 

0.619 
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